Aprendizado de Máquina I

Hugo Tremonte de Carvalho

hugo@dme.ufrj.br


03 - Validação Cruzada

Relembrando um pouco…

Formulação matemática

  • \(Y \in \mathbb{R}\) variável de resposta

  • \(\mathbf{X} = (X_1, \dots, X_p) \in \mathbb{R}^p\) conjunto de \(p\) preditores para \(Y\)

  • Forma geral da relação entre \(\mathbf{X}\) e \(Y\): \[Y = r(\mathbf{X}) + \varepsilon\]

    • \(r\) função fixa porém desconhecida
    • Informação sistemática que \(\mathbf{X}\) dá sobre \(Y\)
    • \(\varepsilon\) erro de observaçao (independente de \(\mathbf{X}\) e de média zero)

Formulação matemática

  • Objetivos: Estimar \(r\) a partir de um conjunto de observações \((\mathbf{x}_i, y_i)_{i = 1, \dots, n}\), bem como aferir a qualidade de tal estimativa

  • Importante: Assumiremos observações i.i.d., ou seja, \((\mathbf{X}_1, Y_1), \dots, (\mathbf{X}_n, Y_n) \sim (\mathbf{X}, Y)\)

    • O que não significa que a distribuição conjunta entre \(\mathbf{X}\) e \(Y\) será conhecida!
    • Em um mundo ideal, temos \(r(\mathbf{x}) = \mathbb{E}[Y | \mathbf{X} = \mathbf{x}]\)
  • \(Y\) qualitativa: problema de regressão

Formulação matemática

  • Objetivo preditivo: Como podemos criar uma função \[g: \mathbb{R}^p \to \mathbb{R}\] que tenha bom poder preditivo? Isto é, como criar \(g\) tal que, dadas novas observações i.i.d. \[(\mathbf{X}_{n + 1}, Y_{n + 1}), \dots, (\mathbf{X}_{n + m}, Y_{n + m}) \sim (\mathbf{X}, Y)\] tenhamos \[g(\mathbf{x}_{n + 1}) \approx y_{n + 1}, \dots, g(\mathbf{x}_{n + m}) \approx y_{n + m}?\]

A função risco

Critério para medir o desempenho de uma dada função de predição \(g: \mathbb{R}^p \to \mathbb{R}\)

  • \((\mathbf{X}, Y)\) nova observação que não foi usada para se encontrar \(g\)
  • \(L(g; (\mathbf{X}, Y))\) – função perda
    • \(L(g; (\mathbf{X}, Y)) = (Y - g(\mathbf{X}))^2\) – perda quadrática
    • \(L(g; (\mathbf{X}, Y)) = |Y - g(\mathbf{X})|\) – perda absoluta
    • muitas outras…

A função risco

  • Tirando a dependência na observação específica…
    • \(R_{\mathrm{pred}}(g) = \mathbb{E}[L(g; (\mathbf{X}, Y))]\)
    • No caso da perda quadrática (usual em regresão) \[R_{\mathrm{pred}}(g) = \mathbb{E}[(Y - g(\mathbf{X}))^2]\]
  • Na prática, estimamos \(R_{\mathrm{pred}}(g)\) através do erro quadrático médio: \[\mathrm{EQM} = \frac{1}{n} \sum_{i = 1}^{n} (y_i - g(\mathbf{x}_i))^2\]

Trazendo a teoria para os dados

Estimando o risco

  • Teorema: Se \[(\mathbf{X}_{n + 1}, Y_{n + 1}), \dots, (\mathbf{X}_{n + m}, Y_{n + m}) \sim (\mathbf{X}, Y)\] são novas observações i.i.d. (ou seja, não foram usadas para se estimar \(g\)), então o EQM calculado em tais observações é um estimador não-viesado para \(R_{\mathrm{pred}}(g)\).

  • Obs.: O valor esperado é tomado somente nas novas observações, e não no conjunto de dados usado para se estimar \(g\).

Estimando o risco

  • Qualitativamente:

    • EQM pequeno \(\iff\) previsões próximas às respostas verdadeiras
    • EQM grande \(\iff\) algumas previsões diferem bastante das respostas verdadeiras
  • Mas onde esse EQM deve ser calculado? Com que conjunto de dados?

  • Obs.: Paralelo entre “aferir qualidade de um modelo” e “selecionar o melhor modelo”

Estimando o risco

  • Conjunto de treinamento: \((\mathbf{x}_{k}, y_{k})_{k = 1, \dots, n}\), usado para estimar a função \(g\)
  • Conjunto de teste: \((\mathbf{x}_{n + k}, y_{n + k})_{k = 1, \dots, m}\), usado para estimar a qualidade de \(g\), através do EQM de teste: \[\frac{1}{m}\sum_{i = n + 1}^{n + m} (y_i - g(\mathbf{x}_i))^2\]
  • Aparentemente queremos \(g\) que minimize o EQM de teste!
  • NÃO é razoável assumir que os EQM de treinamento e teste são semelhantes!

Um exemplo

Figura 2.9 de [ITSL]

Esquerda: Dados simulados de \(r\) (em preto), e três estimativas \(g\) com diferentes complexidades. Direita: EQM de treinamento (cinza) e de teste (vermelho), em função da complexidade do modelo.

Um exemplo mais concreto

  • \(X \sim \mathcal{U}[0, 2]\)
  • \(f(x) = \sin(\pi x)\)
  • \(Y = f(X) + \varepsilon\)
  • \(\varepsilon \sim \mathcal{N}(0, 0.4)\)
  • Ajuste polinomial para diversos valores de \(p\): \[g(x) = \beta_0 + \sum_{n = 1}^{p} \beta_n x^n\]
  • Fonte: Bias-variance trade-off with Python example

Dados gerados segundo o modelo

Ajustes polinomiais

Best fit polynomial degree:  7

Voltando à estimativa do risco

  • Classe de modelos “candidatos”: \(\mathbb{G}\)
  • Modelo \(g\) com baixo EQM de treinamento não é um bom sinal de ajuste aos dados, mas também ao ruído subjacente às observações - sobreajuste ou overfitting
  • Modelo \(g\) com menor EQM de teste também não é o ideal: super-ajuste ao conjunto de teste! Isto é, isso pode levar à escolha de uma função \(g\) que tenha bom desempenho preditivo no conjunto de teste, mas não em novas observações - mesmas observações usadas para comparar todas os \(g\)!

Problema - Como estimar bem o risco e selecionar o melhor modelo?

  • Conjunto de dados \((\mathbf{x}_{k}, y_{k})_{k = 1, \dots}\), i.i.d. segundo \((\mathbf{X}, Y)\)
  • \(\mathbb{G}\) classe de modelos “candidatos” para ajustar aos dados
  • Aferir “bem” a capacidade preditiva para “toda” \(g \in \mathbb{G}\)
  • Escolher a \(g \in \mathbb{G}\) com melhor capacidade de generalização fora do conjunto de dados inicial
  • Escolher o melhor hiperparâmetro para alguma \(g\)
    • Exemplos de hiperparâmetros: \(\lambda\) no Lasso ou Ridge, grau \(p\) de um ajuste polinomial, etc…

Validação cruzada

Dividindo bem o conjunto de dados…

\[\underbrace{(\mathbf{x}_1, y_1), \dots, (\mathbf{x}_n, y_n)}_{\text{conjunto de treinamento}}, \underbrace{(\mathbf{x}_{n + 1}, y_{n + 1}), \dots, (\mathbf{x}_{n + m}, y_{n + m})}_{\text{conjunto de teste}}\]

  • Aproximadamente 70% para treinamento e 30% para teste
    • Treinar modelo requer mais informação \(\Rightarrow\) precisa de mais dados
    • Para ter uma boa estimativa do risco no conjunto de teste, menos dados são necessários

…e fazendo bom uso deles!¹

¹ Especialmente do conjunto de treinamento!

Figura 5.5 de [ITSL]

Ilustração esqumática da validação cruzada com \(k\) dobras

Validação cruzada com \(k\) dobras

  • Randomizar o conjunto de teste (por quê?)
  • Dividir em \(k\) partes iguais (“dobras”, “fold”)
  • Primeira dobra: conjunto temporário de teste
  • \(k - 1\) outras dobras: conjunto temporário de treinamento
    • Treina um modelo \(g_1\) e calcula seu EQM na primeira dobra
  • Repete: \(k\) observações de EQM

Validação cruzada com \(k\) dobras

  • \(\mathrm{EQM}_i\), para \(i = 1, \dots, k\)
    • Modelos \(g_i\) treinados no azul e testados no bege, para \(i = 1, \dots, k\)
  • \(\mathrm{EQM}_{k-\mathrm{VC}} = \frac{1}{k} \sum_{i = 1}^{k} \mathrm{EQM}_i\)
  • Atenção: Não mede a qualidade de um ajuste \(g\) mas sim de um modelo \(g\). Qual a diferença?!

Como escolher o valor de \(k\)?

  • \(k = 1\): treina no conjunto de treinamento e aceita
    • “Treina no treino e testa no teste”: super-ajuste ao conjunto de teste
  • \(k = n\): Leave One Out Cross Validation (LOOCV)
    • Requer \(n\) ajustes do modelo; computacionalmente caro
  • Tipicamente: \(k\) entre 5 e 10
    • Por quê?
    • A resposta usa o balanço entre o viés e a variância

Como escolher o valor de \(k\)?

  • LOOCV
    • Treina \(n\) modelos, cada um em \(n - 1\) amostras
    • Modelos mais “bem treinados” \(\Rightarrow\) mais próximos da “verdade” \(\Rightarrow\) menor viés
    • Conjuntos com grande interseção entre si \(\Rightarrow\) modelos muito correlacionados \(\Rightarrow\) \(\mathrm{EQM}_{n-\mathrm{VC}}\) com alta variância
    • O valor de \(\mathrm{EQM}_{n-\mathrm{VC}}\) é menos confiável!

Como escolher o valor de \(k\)?

  • \(k\)-VC, com \(k \ll n\)
    • Conjuntos com menos interseção entre si \(\Rightarrow\) modelos menos correlacionados \(\Rightarrow\) \(\mathrm{EQM}_{k-\mathrm{VC}}\) com menor variância \(\Rightarrow\) valor mais confiável
    • E quanto ao viés? Se o conjunto de treinamento é “suficientemente grande” e o modelo não é “muito complexo”, em cada etapa o modelo provavelmente estará bem treinado \(\Rightarrow\) baixo impacto no viés!
  • Portanto, \(k\) entre 5 e 10 faz sentido!

E para escolher o melhor ajuste?

\[\underbrace{(\mathbf{x}_1, y_1), \dots, (\mathbf{x}_n, y_n)}_{\text{conjunto de treinamento}}, \underbrace{(\mathbf{x}_{n + 1}, y_{n + 1}), \dots, (\mathbf{x}_{n + m}, y_{n + m})}_{\text{conjunto de teste}}\]

  • Escolhe(m)-se bom(ns) candidato(s) \(g\) conforme o procedimento anterior
  • Treina(m)-se tal(is) candidato(s) \(g\) em todo o conjunto de treinamento. Por quê?
  • Para o(s) modelo(s) \(g\) conforme acima, testa-se o(s) seu(s) desempenho(s) no conjunto de teste
  • Escolhe-se o melhor ajuste com base em tal desempenho

Continuando no exemplo anterior

Relembrando o exemplo

  • \(X \sim \mathcal{U}[0, 2]\)
  • \(f(x) = \sin(\pi x)\)
  • \(Y = f(X) + \varepsilon\)
  • \(\varepsilon \sim \mathcal{N}(0, 0.4)\)
  • Como escolher o melhor valor de \(p\) no ajuste polinomial? \[g(x) = \beta_0 + \sum_{n = 1}^{p} \beta_n x^n\]

Implementando o procedimento

  • Gerando dados já separados por treinamento e teste

    n_obs_tr = 300
    n_obs_vl = 100
    degs = [0, 3, 5, 7, 10, 15, 20]
    
    x_tr, y_tr = generate_data(n_obs = n_obs_tr, err_stdv = np.sqrt(0.4))
    x_vl, y_vl = generate_data(n_obs = n_obs_vl, err_stdv = np.sqrt(0.4))
  • Separando o conjunto de treinamento em \(k\) dobras

    kf = KFold(n_splits = 5, shuffle = True)
    for i, (train_index, test_index) in enumerate(kf.split(x_tr, y_tr)):
      print(f"Fold {i}:")
      print(f"  Train: index={train_index}")
      print(f"  Test:  index={test_index}")
    Fold 0:
      Train: index=[  1   2   4   5   6   7   8   9  12  13  14  16  17  18  20  21  23  24
      25  26  27  28  29  31  32  34  35  36  37  38  39  40  42  43  46  47
      48  49  50  51  52  53  54  55  56  59  60  61  62  63  64  66  68  69
      70  71  72  73  74  75  77  78  79  80  81  82  83  84  85  86  87  88
      89  90  92  93  94  95  96  99 101 102 103 104 105 107 108 109 110 112
     113 114 115 116 117 118 119 121 122 125 127 128 129 130 131 132 133 135
     136 137 138 139 140 141 143 145 146 149 150 151 152 153 154 155 156 157
     158 159 160 161 162 163 164 166 167 168 169 170 171 173 175 176 178 180
     181 184 185 186 187 188 189 190 193 194 195 196 197 198 199 201 204 205
     206 208 210 211 212 213 214 215 216 217 218 219 221 222 223 225 226 227
     228 229 230 231 232 233 235 236 237 238 239 240 241 242 244 245 247 248
     250 251 252 253 255 257 258 259 260 261 262 263 265 266 268 269 270 271
     272 273 274 276 277 278 279 281 282 283 285 286 287 288 289 290 291 292
     293 294 295 296 298 299]
      Test:  index=[  0   3  10  11  15  19  22  30  33  41  44  45  57  58  65  67  76  91
      97  98 100 106 111 120 123 124 126 134 142 144 147 148 165 172 174 177
     179 182 183 191 192 200 202 203 207 209 220 224 234 243 246 249 254 256
     264 267 275 280 284 297]
    Fold 1:
      Train: index=[  0   1   2   3   4   5   6   7   8   9  10  11  12  14  15  17  18  19
      21  22  23  24  25  27  28  29  30  31  32  33  35  37  38  39  40  41
      42  43  44  45  46  47  48  50  51  52  53  54  55  56  57  58  59  60
      61  62  63  65  66  67  71  74  76  77  78  80  81  82  84  85  86  88
      90  91  94  95  96  97  98  99 100 101 102 103 105 106 107 108 109 110
     111 113 114 115 117 119 120 121 122 123 124 125 126 127 129 130 131 132
     133 134 135 136 141 142 143 144 145 146 147 148 149 150 152 153 154 155
     156 158 159 161 162 163 164 165 166 167 168 169 171 172 174 176 177 178
     179 180 181 182 183 184 185 186 188 189 190 191 192 193 194 195 196 197
     198 200 202 203 204 206 207 208 209 211 213 214 216 219 220 222 223 224
     225 226 227 228 229 230 231 234 235 236 237 238 239 240 242 243 245 246
     248 249 251 252 253 254 255 256 258 259 260 261 262 264 265 267 268 269
     270 271 272 273 274 275 277 278 279 280 281 282 283 284 285 286 287 288
     289 290 292 296 297 298]
      Test:  index=[ 13  16  20  26  34  36  49  64  68  69  70  72  73  75  79  83  87  89
      92  93 104 112 116 118 128 137 138 139 140 151 157 160 170 173 175 187
     199 201 205 210 212 215 217 218 221 232 233 241 244 247 250 257 263 266
     276 291 293 294 295 299]
    Fold 2:
      Train: index=[  0   1   2   3   4   6   8  10  11  12  13  14  15  16  17  18  19  20
      22  23  24  26  27  28  29  30  31  32  33  34  35  36  37  40  41  42
      43  44  45  47  49  52  53  57  58  59  60  61  62  64  65  66  67  68
      69  70  72  73  74  75  76  77  79  80  81  82  83  84  85  86  87  88
      89  90  91  92  93  94  95  97  98  99 100 101 102 103 104 106 108 109
     110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 128
     129 130 133 134 135 137 138 139 140 141 142 143 144 145 146 147 148 149
     150 151 152 154 156 157 159 160 161 162 163 165 167 168 169 170 172 173
     174 175 176 177 178 179 180 182 183 185 187 189 190 191 192 193 194 195
     198 199 200 201 202 203 204 205 206 207 209 210 211 212 214 215 216 217
     218 219 220 221 222 223 224 225 228 231 232 233 234 237 238 241 243 244
     245 246 247 248 249 250 251 252 253 254 255 256 257 258 260 262 263 264
     265 266 267 269 270 272 274 275 276 279 280 281 282 283 284 286 290 291
     292 293 294 295 297 299]
      Test:  index=[  5   7   9  21  25  38  39  46  48  50  51  54  55  56  63  71  78  96
     105 107 127 131 132 136 153 155 158 164 166 171 181 184 186 188 196 197
     208 213 226 227 229 230 235 236 239 240 242 259 261 268 271 273 277 278
     285 287 288 289 296 298]
    Fold 3:
      Train: index=[  0   2   3   5   7   9  10  11  13  14  15  16  17  18  19  20  21  22
      23  24  25  26  27  28  29  30  31  32  33  34  36  37  38  39  40  41
      43  44  45  46  48  49  50  51  54  55  56  57  58  59  63  64  65  66
      67  68  69  70  71  72  73  74  75  76  77  78  79  82  83  84  85  87
      89  90  91  92  93  94  96  97  98  99 100 104 105 106 107 108 109 110
     111 112 113 114 115 116 117 118 119 120 121 123 124 126 127 128 131 132
     134 136 137 138 139 140 141 142 144 146 147 148 151 152 153 154 155 156
     157 158 159 160 162 164 165 166 170 171 172 173 174 175 177 179 180 181
     182 183 184 186 187 188 189 191 192 194 196 197 198 199 200 201 202 203
     204 205 206 207 208 209 210 212 213 214 215 217 218 219 220 221 223 224
     226 227 228 229 230 231 232 233 234 235 236 237 239 240 241 242 243 244
     246 247 249 250 254 256 257 258 259 260 261 262 263 264 265 266 267 268
     271 273 274 275 276 277 278 280 284 285 286 287 288 289 290 291 292 293
     294 295 296 297 298 299]
      Test:  index=[  1   4   6   8  12  35  42  47  52  53  60  61  62  80  81  86  88  95
     101 102 103 122 125 129 130 133 135 143 145 149 150 161 163 167 168 169
     176 178 185 190 193 195 211 216 222 225 238 245 248 251 252 253 255 269
     270 272 279 281 282 283]
    Fold 4:
      Train: index=[  0   1   3   4   5   6   7   8   9  10  11  12  13  15  16  19  20  21
      22  25  26  30  33  34  35  36  38  39  41  42  44  45  46  47  48  49
      50  51  52  53  54  55  56  57  58  60  61  62  63  64  65  67  68  69
      70  71  72  73  75  76  78  79  80  81  83  86  87  88  89  91  92  93
      95  96  97  98 100 101 102 103 104 105 106 107 111 112 116 118 120 122
     123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140
     142 143 144 145 147 148 149 150 151 153 155 157 158 160 161 163 164 165
     166 167 168 169 170 171 172 173 174 175 176 177 178 179 181 182 183 184
     185 186 187 188 190 191 192 193 195 196 197 199 200 201 202 203 205 207
     208 209 210 211 212 213 215 216 217 218 220 221 222 224 225 226 227 229
     230 232 233 234 235 236 238 239 240 241 242 243 244 245 246 247 248 249
     250 251 252 253 254 255 256 257 259 261 263 264 266 267 268 269 270 271
     272 273 275 276 277 278 279 280 281 282 283 284 285 287 288 289 291 293
     294 295 296 297 298 299]
      Test:  index=[  2  14  17  18  23  24  27  28  29  31  32  37  40  43  59  66  74  77
      82  84  85  90  94  99 108 109 110 113 114 115 117 119 121 141 146 152
     154 156 159 162 180 189 194 198 204 206 214 219 223 228 231 237 258 260
     262 265 274 286 290 292]

Implementando o procedimento

mse_deg = {}

for i, (train_index, test_index) in enumerate(kf.split(x_tr, y_tr)):
  x_tr_aux = x_tr[train_index]
  y_tr_aux = y_tr[train_index]
  x_te_aux = x_tr[test_index]
  y_te_aux = y_tr[test_index]

  fit_coeffs = fit_poly_models(degs, x_tr_aux, y_tr_aux)
  y_te_aux_est = eval_poly_models(fit_coeffs, x_te_aux)
  mse_te_aux = comp_poly_mse(y_te_aux_est, y_te_aux)
  for k, v in mse_te_aux.items():
    if k in mse_deg:
      mse_deg[k].append(v)
    else:
      mse_deg[k] = [v]
print(mse_deg)
{'0': [1.0835319737812303, 0.9123508710126779, 0.8515078166111659, 0.8211992644333059, 0.906370276956641], '3': [0.4685550190713497, 0.4413105599865085, 0.5172188314619182, 0.38201835586616856, 0.3723633515396271], '5': [0.4463387593756728, 0.43852306378300354, 0.5256267862112263, 0.3826292538478435, 0.36739045422128197], '7': [0.4597438072888968, 0.43815816190141615, 0.5314811538962878, 0.3994972398238271, 0.3752576591112432], '10': [0.4586051026196262, 0.44153075197022534, 0.5332320817922022, 0.412506993610928, 0.40613345148146784], '15': [0.4599156163082655, 0.4551540655533446, 0.5472345688205718, 0.4140035858455774, 0.42128232100939694], '20': [0.45284758279986315, 0.46718928191951736, 0.586792317822178, 0.6681552868787054, 0.40650534309922826]}

Implementando o procedimento

mse_deg_mean = {}
for k, v in mse_deg.items():
  mse_deg_mean[k] = np.mean(np.array(v))

print(mse_deg_mean)

print("Best fit polynomial degree (OLD): ", degs[list(mse_vl.values()).index(min(list(mse_vl.values())))])

print("Best fit polynomial degree (NEW): ", degs[list(mse_deg_mean.values()).index(min(list(mse_deg_mean.values())))])
{'0': 0.9149920405590042, '3': 0.43629322358511435, '5': 0.4321016634878056, '7': 0.4408276044043342, '10': 0.45040167629488986, '15': 0.45951803150743126, '20': 0.5162979625038984}
Best fit polynomial degree (OLD):  7
Best fit polynomial degree (NEW):  5